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We study the adsorption of primitive model electrolytes into a layered slit system using grand canonical Monte 
Carlo simulations. The slit system contains a series of charged membranes. The ions are forbidden from the 
membranes, while they are allowed to be adsorbed into the slits between the membranes. We focus on the 
electrical properties of the slit system. We show concentration, charge, electric field, and electrical potential 
profiles. We show that the potential difference between the slit system and the bulk phase is mainly due to the 
double layers formed at the boundaries of the slit system, but polarization of external slits also contributes to 
the potential drop. We demonstrate that the electrical work necessary to bring an ion into the slit system can 
be studied only if we simulate the slit together with the bulk phases in one single simulation cell. 
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1. Introduction 

Electrical double layers (DLs) formed in an electrol5fte near a charged surface have been studied using 
both simulation |1-12] and theoretical |13- 23] methods. The basis of these studies is the Primitive Model 
(PM) of electrolytes that represents the solvent as a dielectric continuum. Simulations necessarily use a 
finite simulation cell for these studies, where the electrolyte is confined between two charged walls. 

If these walls are far enough from each other, the DLs formed at the walls are independent of each 
other and a charge neutral bulk region is formed in the middle of the simulation cell. The reference point 
for the electrical potential then can be set in this bulk region. If the walls, however, are close to each 
other so that the two DLs overlap (slit), the bulk region disappears |24- 34]. Slits are generally simulated 
in the grand canonical (GC) ensemble, where the electrolyte in the slit is in equilibrium with a virtual 
bulk phase represented by its temperature and the chemical potentials of the ionic species. The ground 
of the electrical potential cannot be set in this bulk phase, because it is not present in the simulation cell 
and Poisson's equation cannot be integrated over it. 

The studies for electrolytes confined in a sht considered only a single slit (we will call it the "lonely 
slit") in equilibrium with a bulk in the GC ensemble [24-34]. In this work we are primarily concerned 
with the electrical properties of the sUt. Most importantly, we want to calculate the potential difference 
between a slit and the bulk phase. Namely, we are concerned with the electrical work needed to bring 
an ion from bulk into the slit. Therefore, we simulate the slits confined between membranes with bulk 
electrolytes on the other sides of the membranes. 

Furthermore, we consider more than one slit (we will call these "slit systems"). Slit systems are first- 
order models for layered silicate minerals |35- 38], porous electrodes |39-43], and lyotropic lamellar liq- 
uid crystals (T^H). The structure, swelling, and adsorption properties of such materials (e.g. kaolinite, 
montmorillonite) are subject of extensive experimental and simulation studies fii'-^ii']. In this work, we 
focus on the electrical properties of charged slit systems and show concentration, charge, electric field. 
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and potential profiles for different geometrical parameters (width of the slit, width of the membrane), 
membrane charge, and electrolytes (concentration, ionic charge). 



2. Model of the slit system 

The electrolyte is modeled with the Restricted Primitive Model (RPM). In this model, the solvent is 
represented by its dielectric response characterized by the dielectric constant e, while the ions are repre- 
sented by charged hard spheres interacting through the screened Coulomb + hard sphere pair potential: 

(oo for rij < Ri + Rj , 

^±31 for ru^R^+Rj, ^^.l) 
4:/reoe- rtj 

where rij is the distance of two ions, qi is the ionic charge (qi - Zj-e, Zi being the valence and e the 
elementary charge), and Rj is the ionic radius. In the RPM, Ri is the same for every ionic species (R - 1.5 
A in this work). The ionic charges are point charges in the centers of hard spheres. It was shown that this 
simple model can reproduce the non-monotonic concentration dependence of the activity coefficient of 
electrolytes as soon as we assume that the dielectric constant is concentration dependent (s^. Since this 
work is a model calculation, we do not change the dielectric constant with the concentration, but we fix 
it at the value e = 78.4. 




Figure 1 . Geometry of the slit system and results for a 0.5 M 1:1 electrolyte in a slit system with parameters 
JVs = 3, Ls = 10 A, Lm = 5 a, and ctm = -0.05 Cm~^. The concentration and charge profiles are measured 
in M, while the electric field profile is measured in . Here and in later figures, the membrane regions 
are indicated with shaded grey areas. 
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The ath membrane is confined by two ttard walls at xjj and x^, where each hard wall can carry a ctm 
surface charge. The interaction potential between such a charged hard wall and an ion is 

{oo for |x| < Ri , 

-"'""^'"1 for \x\^Ri, (2.2) 

where \x\ is the distance of the ion from the surface. 

There are Nm membranes of width Lm (this is the distance of the two walls forming the membrane: 
Lm - for every a) in the slit system . In this work, both membrane walls carry a ctm surface charge. 

The distance of two membranes, namely, the width of the slit is Ls- This distance is kept fixed during the 
simulation, namely, the slit system is rigid. We should keep in mind, however, that the overlapping DLs 
between the like-charged macromolecules play a central role in the attractive force acting between them 
and contribute to the cohesion of such materials 138, 31- 531. 

Ions are not allowed to enter the membranes, so the following Boltzmann factor is used to forbid the 
ions from the membranes: 



exp 



kT 



for x^ + Ri<x<xl-Ri, 3) 

1 otherwise, 



where x}^ and x^ are the coordinates of the left and right walls of the ath membrane. Obviously, the 
number of slits is - Nm - 1. 

There are two bulk phases of widths Lb on the two sides of the slit system. The simulation cell is closed 
by two hard walls on the left hand side of the left bulk region and on the right hand side of the right bulk 
region. In this study, these walls are uncharged. The geometry is illustrated in figure[T] 

3. Grand Canonical Monte Carlo simulations 

Grand Canonical Monte Carlo (GCMC) simulations have been performed for the system described 
above. Periodic boundary conditions have been applied in the y - z dimensions. The effect of the ions 
outside the central simulation cell has been taken into account by the charged sheet method proposed by 
Torrie and Valleau |2] and developed further by Boda et.al. |54]. 

In GCMC simulations of the DL, in addition to the usual particle displacement steps, we insert and 
delete neutral clusters of ions, e. g., v+ cations and v_ anions (v+ and v_ being the stoichiometric coeffi- 
cients). This way, we make sure that the simulation cell is charge neutral in every instant of the simula- 
tion. The acceptance probabUity of these steps is 

where N+ and N- are the numbers of cations and anions before insertion/deletion, V is the volume into 
which we insert the centres of the ions, v - v+ + v-, x = 1 for insertion and x = -1 for deletion, while 
At/ is the energy change of the insertion/deletion that contains the pair energies between ions, ions and 
charged walls, and ions and charged sheets (it becomes infinite in the case of overlap). The quantity 
fi± - v+/i+ -I- v_/i_ is the chemical potential of the salt (/x+ and /i_ being the chemical potentials of the 
ions) determined by the Adaptive-GCMC (A-GCMC) method of Malasics et al. (isjliill. 

The dimensions of the simulation cell in the y-z dimensions is in the range of 120-150 A. System-size 
checks indicated little sensitivity of the concentration profiles on the y-z dimensions of the cell. In a 
typical simulation, the sample contained several hundreds of mUUons (10^) configurations. 

4. Solution of Poisson's equation 

The main output quantities of the simulations are the density profiles of the various ionic species. 
Pi (x), from which the ionic charge profile is obtained as 

^(x) = ^Ziep,W. (4.1) 
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Poisson's equation 



1 



-qtotix) 



(4.2) 



is solved for the mean electrical potential <l)(x) using Neumann boundary conditions. The total charge 
density qtot ix) contains all charges in the system (including the ctm membrane charges in addition to the 
ionic charge density qixj): 

qtot{x)^q{x) + Y,(JiS{x-Xi], (4.3) 



where at surface charges are placed in positions x/. After integration we obtain 

X X 

dO(x) 



dx 



-— f qtotix')dx' + Ci^ — — f q[x')dx' - — Y Ui + Ci, 
eoe J eoe J 



(4.4) 



where the notation /(x, < x) under the sum means that we include only those sheets in the sum whose 
coordinates are smaller than x. The numerical integration is performed using the rectangle rule. The 
electric field is defined as 

£(x) = r-^. (4.5) 



dx 



Since our system is charge neutral: 



oo 

j t7(x)dx + ^(TiA = 0, 



(4.6) 



the electric field outside the simulation domain is zero {A - is the area of the simulation cell in the 
y- z plane). The first boundary condition, therefore, is 



£(-00) = 0, 



(4.7) 



from which the value Ci = is obtained for the first integration constant. After one more integration the 
electrical potential is obtained 



<D(x) = 



X X 

— [ [ ^tot(:t")dx" 
eoe J J 

—00 L— 00 

X T x' 

— [ [ q[x")dx" 
eoe J J 



dx' + C1X + C2 



L '^i 

^oe i(xi<x) 



C2, 



(4.8) 



where we assumed that Ci - 0. Numerical integration is performed according to the trapezoidal rule. 
Since the system is charge neutral, the Neumann boundary condition on the other side (£(00) = 0) is 
automatically fulfilled. 

The zero level of the potential is arbitrary. Therefore, we choose the integration constant C2 so that 
the potential is zero in the left hand side bulk. This is done by starting with C2 - 0, averaging the potential 
over the left hand side bulk (where it is constant), and deducting this value from the potential profile. 

Note that this method of solving the Poisson's equation is different from the traditional convolution 
integral 



0(x) = / (x'-x)f7tot(x')dx', 

eoe J 



(4.9) 



where the boundary condition is set in a bulk region (Loo). This formula is used in theoretical studies of 
an isolated DL, where the boundary condition is set infinitely far from the electrode (Loo ^ 00, where 
the potential and the electric field are zero). This equation can be used in simulation studies if there is a 
bulk region in the simulation cell and if the statistics is good enough. If there is no bulk region, such as 
in the case of a lonely slit, the convolution integral is not practical. Furthermore, equation <4.8> provides 
results for the potential profiles with much better accuracy than the convolution integral. A formula that 
is equivalent to equation i4.8\ was introduced by Kiyohara and Asaka 13 9.1 without any reference to the 
fact that it corresponds to the Neumann boundary conditions. 
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5. Results and discussion 

In this work, we measure the distances in A, so particle densities are measured in A~^. In the figures, 
however, we plot concentration profiles that are related to the density profiles through c, (x) - 1660.58 x 
Piix) (the unit of concentration is mol/dm^ abbreviated as M). The charge profUe is also computed in 
terms of concentrations and is normaUzed by the elementary charge: q{x) le - z+ c+ (x) + z_ c_ (x) (the unit 
is M). The potential is plotted in a dimensionless form as eO(x)/fcr. The electric field is the derivative of 
this dimensionless potential, so its unit is A~^. We wiU denote it by E* (x). 
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Figure 2. (Color online) The effect of the number of slits for a 0.5 M 1:1 electrolyte in slit systems with 
parameters Lg = 10 A, Ljvi = 5 A, and = -0.05 Cm~^. Results for = I and 3 are shown. For the 
Ws = 1 case, results are also shown for the lonely slit, when the bulk regions are not included in the 
simulations (symbols). The outermost membrane region (shaded with Ughter grey) is present only for the 
case ATg = 3, so they are not forbidden regions for the solid red line concentration profiles for Ns = 1 that 
seem to overlap with them. The role of the symbols is to distinguish the curve from other, overlapping 
curves; they indicate selected simulation points. This is also true for other figures, where symbols are 
used. 

The temperature is T = 298.15 K. We used 1:1 and 2:1 RPM electrolytes with ions of radii R - 1.5 A. 
The size of the bulk phases was Lb = 50 A (simulations with Lb = 80 A gave the same results). We show 
the results for various values of the number of slits (TVs), width of slits (Ls), width of membranes (Lm), 
membrane charge (ctm). and electrolyte concentrations (c). 

A typical result is seen in figure [l] where concentration, charge, electric field, and potential profUes 
are plotted for a 1:1 electrolyte adsorbed in a slit system with parameters indicated in the caption. The 
concentration profUes of the cations (counter-ions) show peaks in the slits that are necessary to balance 
the membrane charges. At these parameters, there is enough space in the slits to accept some colons 
too. The concentration profiles of both ions approach the bulk values in the two bulk regions on the two 
sides of the slit system. Near the external membranes double layers are formed. These double layers are 
apparent in the charge profile. They also balance the membrane charge but they extend into the bulk, 
because space is available to form diffuse layers. 
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These double layers are responsible for the poten- 
tial drops between the bulk regions and the outmost 
membranes of the slit system (Odl. see figure [!)• Let 
us define the membrane potential, as the potential 
difference between the inner membranes and the bulk. 
The potential drop associated with the external DLs, 
Odl. dominates the membrane potential, <I>m- The dif- 
ference is due to the fact that there is less ionic charge 
in the outer shts than necessary to balance the mem- 
brane charge. The missing charge in the outer DLs bal- 
ances the membrane charges from outside. 

The fact that the outer slits and the associated mem- 
brane charges do not cancel is proved by the electric Figure 3. The counter-ion concentration pro- 
field profile. The electric field is zero in the inner mem- files of figure|2]in the innermost slit, 
branes because inside the slit system the adsorbed ionic 

charge totally balances the membrane charge. In the outer membranes, however, the electric field is non- 
zero indicating that there are more counter-ions in the outer DLs and less counter-ions in the outer slits. 
(The non-zero electric field is reflected by the slope of the potential in the outer membranes.) This polar- 
ization adds an additional term to Om denoted by <1>pol- Of course, the relation 3>dl + 'I'pol = "I'm between 
these potential values holds. In this case, OpoL is much smaller than Odl- 

The effect of the number of slits Figure [2] shows the results for different values of Ns (Afs - 1 and 3), 
the lonely slit (A/s - 1, symbols), among them. Only the outermost slit is polarized: the inner slits are not 
effected by the outer world. The concentration and the potential profiles show a periodic behaviour in 
the inner slits. This fact is even better demonstrated in later figures where we show the results for more 
slits. 

An interesting result is that the behaviour of the lonely slit is quite similar to the behaviour of the in- 
nermost slit when bulk phases are present. The shape of the potential profile is similar, but the reference 
point of the potential profile is not well defined. At least, this cannot be related to the bulk region because 
the bulk region is not present in the simulation cell when the lonely slit is simulated. 

The shape of concentration profiles also look similar in the lonely slit and in the innermost slits, at 
least, at the scale of figure[2] Figure [3] magnifies the peaks of the counter-ion profiles. The ionic profiles 
in the lonely slit show a maximum at contact position: a typical hard sphere effect at the wall. The ionic 
profiles in the slit system, on the other hand, show a smooth decrease at the wall. This is an electrostatic 
effect: the dipole field of the polarized charges outside the slit system exerts a repulsive effect on the 
counter-ions. 

This result shows that the lonely sUt is a good approximation if we want to study the structure of the 
density profiles in the slit, but this is insufficient if we want to get information about the electrical work 
that is needed to bring an ion into the sht. In that case, the sht should be simulated together with the bulk 
phases outside it. 

The effect of the width of slits In the slits, the potential profile declines (in absolute value) to zero, but 
it does not reach zero because charge neutral bulk regions are not formed in the narrow slits. If the slits 
were wide enough, independent DLs would form near the membranes and the potential would drop to 
zero. Figure|4]shows this effect with an increasing slit width. 

As Ls increases, the average concentration in the slits decreases and the cation and anion profiles in 
the middle of the slits become more and more similar. The potential profile, therefore, gets closer to zero. 
The potential drop in the DL, <I>dl. is the same in aU cases: the interior of the slit system does not effect the 
outer DL. The potential profile inside the slit system, on the other hand, sensitively depends on Ls. For 
large 1$, the <1>pol potential drop vanishes, and the potential is the same in all membranes. Polarization 
of the outer membranes (Opol) appears for narrow slits. 
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Figure 5. The effect of the membrane charge for a 
0.5 M 1:1 electroljrte in a sht system with param- 
eters ATg = 6, Ls = 5 A, and = 5 A. Potential 
profiles for values = -0.01, -0.03, -0.05, and 
-0.07 Cm~^ are shown. 



Figure 6. The potential drops between the edge of 
the slit system and the bulk (Odl) and between the 
inner membranes of the slit system and the bulk 
(Om) hased on the potential profiles of figure[5] 
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The effect of the membrane charge Increasing the charge of the membranes, ctm, increases the amount 
of ionic charge that is necessary to balance it. This increases the charge difference between the slits and 
the membranes, which makes the potential fluctuation in the slit system larger (figure [5}. The amount 
of charge in the diffuse layers outside the slit systems also increases, which increases the potential drop 
across these DLs (<1>dl)- The potential increases non-linearly with ctm as seen in figure[6] This phenomenon 
was already observed in the case of DLs 111(1 . 

The DL and POL components of the potential drop are also shown in figure [6] The relative weight of 
the DL component in <1>m is smaller at small membrane charges, while OpoL is relatively small at large 
membrane charges. As a matter of fact, <1>pol shows a non-monotonous behaviour and decreases with 
increasing ctm at large membrane charges. 



The effect of concentration The electrode potential of the DL layer for a given electrode charge (the 
capacitance of the DL) depends on the concentration of the electrolyte. The electrical properties of the 
slit system also depend on the concentration. Figure [7] shows the difference between the results for O.OS 
and 0.5 M electrolytes. The net amount of ionic charge in a slit is fixed by the membrane charge, but this 
is achieved in different ways for the two concentrations. The chemical potential n± that determines how 
many ions tend to be adsorbed in the slit are different at different concentrations. At a higher concen- 
tration (0.5 M), there are more counter-ions in the slit and they drag some colons with themselves. At a 
lower concentration (0.05 M), only counter-ions are in the slit. 




X /A 

Figure 7. The effect of electrolyte concentration for 1:1 electrolytes in a slit system with parameters A/g = 
7, is = 10 A, Im = 10 A, and ctm = -0.02 Cm~^. Concentration, electric field, and potential profiles for 
concentrations c = 0.05 and 0.5 are shown. Only the outer slits at the right hand side of the slit system are 
shown. Anions are indicated with symbols. 



The behaviour of the outer slits is also different in the two cases. There are more missing charges in 
the outermost slit in the 0.05 M case. This is also shown in the non-zero electrical field in the outermost 
membrane. Therefore, the <I>pol potential term is larger at small concentrations. 
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The structure of the outer DLs is also different. As in the case of a DL at an electrode, the $dl potential 
drop is larger at smaU concentrations. 

The effect of divalent ions The overall situation for 2:1 electrolytes is very similar to that for 1:1 elec- 
trolytes. The most notable difference occurs between membrane charges of the opposite sign but the 
same absolute value. Figure |8] shows results for a 2:1 electrolyte in a slit system with ctm = ±0.05 Cm~^ 
membrane charges. 
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Figure 8. The effect of the sign of the membrane charge for a 2:1 electrolyte in a slit system with pa- 
rameters Ns = 10, 1$ = 5 A, and Lm = 5 A. Concentration and potential profiles for membrane charge 
(Tm = ±0.05 Cm~^ are shown. Only the outer slits at the right hand side of the slit system are shown. 
Colons are indicated with symbols. 

The divalent ions are more efficient in balancing the membrane charge than monovalent ions because 
they provide twice the charge occupying approximately the same space. Polarization of the outer sUt, 
therefore, is negligible at a negative membrane charge when the divalent cation is the counter-ion. Since 
the divalent ions form a more compact diffuse layer, the potential drop Odl is smaller in this case. 

At positive membrane charge, on the other hand, both Odl and Opol are larger (in absolute value) 
than in the case of negative membrane charge. 

6. Conclusions 

GCMC simulation results for PM electrolytes adsorbed in a slit system are reported. We conclude that 
it is necessary to study the slit together with the bulk phase with which it is in equilibrium if we want to 
get information on the potential difference between the slit and the bulk. 

This potential difference has two main components. The double layers outside the slit system produce 
potential drops analogous to the electrode potential in the case of a separate DL. This component of the 
potential drop (Odl) depends primarily on the properties of the bulk electrolyte. 

If the ions have difficulty in entering the slits, the outermost slits have less charge than necessary to 
balance the membrane charge. The membrane charge then is neutralized from outside, from the outer 
DLs. This charge missing from the outer shts results in a polarization of these slits and in a potential 
difference. This component of the potential drop (Opol) depends on both the electrolyte properties and 
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geometrical parameters of the slit system. The tpoL componem is larger in absolute value (compared 
to the Odl component) (1) if the width of the slits is smaller, (2) if the electrode charge is smaller in 
absolute value (but not too large), (3) if the concentration is smaller, and (4) if the monovalent ion is the 
counter-ion (the membrane charge positive for 2:1 electrolytes). 

We plan to study (1) competition between counter-ions of different charge and/or size (ionic selectivity 
- this is inspired by our ion channel studies (57|]), (2) ions capable of being adsorbed into the membrane 



regulated by a Boltzmann-factor as in our previous membrane studies 158], and (3) ion transport using 
our Local Equihbrium Monte Carlo method [59,1 . with which a chemical potential gradient can be applied 
across the simulation cell. 
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MoHTe Kapyio ci/iMyyinL^n e/ieKTpi/iHHi/ix BiiacTi/iBOCTeM 
eyieKTpo/iiTiB, aACop6oBaHi/ixy 3ap5iA>KeHi ii4i/ii/iHonoAi6Hi 

CMCTeMI/l 

p. KoBan, M. BayiicKO, fl,. Bofla 

cPaKy/ibTeT c|)i3MHH0T xImIT, VHiBepcMTeT riaHHOHii, Becnpew, Vropmi/iHa 

BuKopucTOByewo MoHie Kap;io CMMy;iHL4ii y se/inKOMy KaHOHiHHOMy aHcawS/ii fl/ifl BUBHeHHJi aACop6L(iy nprnuii- 
TMBHoT MOfle/ii e/ieKipo/iirlB b luapyBaTy inmMHO-nofli6Hy cucreMy. OcraHHA MicrnTb HUSKy sap^flxeHnx mbm- 
6paH. MeM6paHH e HeflOCTynHWMH fl/ia IohIb, ane Moxyib aflCop6yBaTWCb b Mi>KMeM6paHHi mi/iwHW. Mm koh- 
L(eHTpyeMO yBary Ha e/ieKipuHHUx B/iacrnBOcrax mmwHHO'i cucreMW. OrpuMaHO KOHi^eHTpai^iK), sap^fl, e/ie- 
KTpuHHe noyie ra npo<|)i/ib e/ieKipuHHOro noTeHi^ia/iy. noKasaHO, L140 noTeHi4ia;ibHa pi3HWL4fl Mix luI/iuhhok) 
ra o6'£MHOK) c|)a3aMW cnpuHHHeHa, b ocHOBHOMy, BWHWKHeHH^iM Ha rpaHi/mflx mi/iwHW noflBiwHwx liiapiB, ane 
BKJiafl BHOcnTb TaKO>K i no;iflpn3aL4ifl BOBHiiiiHix LL(i;inH. Mvt fleMOHCrpyEMO, mo e;ieKTpw4Ha po6oTa, ^Ky norpi- 
6ho BfliPiCHWTH 3 MeroK) BHeceHHfl iOHa b mi/iwHHy cwcreMy, Moxe 6yTW orpuMana /iwiue iop,\, Konv\ cuMy/iai^ifl 

mi/IMHl/l BUKOHyeibCfl pa30M i3 06'eMHHMW (()a3aMH B OflhiPi CWMy;i51L4iRHiPi K0MipL4i. 

KjiiOMOBi c;iOBa: Mohtb Kapjio, npuMirnBHa MO^enb eneKTponniB, mini/iHi/i 
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